Climate Indices Notebook

This document reproduces the different steps followed before generating the shiny app.

Including relevant libraries

library(tidyverse)  # includes ggplot2 and readr
library(plotly)     # interactive plots

Loading indices file

Mset <- read_delim(file="http://www3.mbari.org/science/upper-ocean-systems/biological-oceanography/GlobalModes/Mset.txt",comment="%",delim="\t")

Parsed with column specification:
cols(
  YEAR = col_integer(),
  MONTH = col_character(),
  M1 = col_double(),
  M2 = col_double(),
  M3 = col_double(),
  M4 = col_double(),
  M5 = col_double(),
  M6 = col_double()
)

Mset$time <- as.Date(paste(Mset$YEAR,Mset$MONTH,rep(15,length(Mset$YEAR))),format="%Y %m %d")
indices <- read_csv(file="http://www3.mbari.org/science/upper-ocean-systems/biological-oceanography/GlobalModes/climate_indices.csv",comment="%")

Parsed with column specification:
cols(
  year = col_integer(),
  month = col_integer(),
  MEI = col_double(),
  PDO = col_double(),
  NPGO = col_double(),
  EMI = col_double(),
  Nino3 = col_double(),
  Nino4 = col_double(),
  AMO = col_double(),
  NAO = col_double(),
  ATL3 = col_double()
)

nb_time <- length(indices$year)
indices$time <- as.Date(paste(indices$year,indices$month,rep(15,nb_time)),format="%Y %m %d")

Plotting MEI & M1

plot(indices$time,indices$MEI,type="l",xlab="Time",ylab="Index")
points(Mset$time,Mset$M1,type="l",col="red")

Testing the use of varname & shading

indexname1 <- "PDO"
indexname2 <- "NPGO"
iok1 <- !is.na(indices[,indexname1])
ipos1 <- indices[,indexname1]>0 & !is.na(indices[,indexname1])
ineg1 <- indices[,indexname1]<0 & !is.na(indices[,indexname1])
x <- indices$time
y0 <- rep(0,nb_time)
plot(x,t(indices[,indexname1]),type="l",xlab="Time",ylab=indexname1)
polygon(c(x[iok1],rev(x[iok1])),c(y0[iok1],rev(t(indices[iok1,indexname1]))),col="skyblue")

points(x,t(indices[,indexname2]),type="l",col="red")

Testing 2 Y-axis

Based on https://www.r-bloggers.com/r-single-plot-with-two-different-y-axes/

indexname1 <- "M1"
indexname2 <- "MEI"
iok1 <- !is.na(Mset[,indexname1])
x <- Mset$time
y0 <- rep(0,length(Mset$YEAR))
par(mar = c(5,5,2,5))
plot(x,t(Mset[,indexname1]),type="l",xlab="Time",ylab=indexname1)
par(new=T)

with(indices,plot(time,MEI,type="l",axes=F,xlab=NA,ylab=NA,col="red",pch=16,cex=1.2))
#plot(indices$time,t(indices[,indexname2]),type="l",axes=F,xlab=NA, ylab=NA,col="red")
axis(side = 4, col="red")

mtext(side = 4, line = 3, indexname2,col="red")
legend("topleft",legend=c(indexname1,indexname2),lty=c(1,1),col=c("black", "red"))

Testing ggplot2 - area plot, double axis, interactive plot

indexname1 <- "M3"
indexname2 <- "PDO"
max1 <- max(Mset[,indexname1],na.rm=TRUE)
max2 <- max(indices[,indexname2],na.rm=TRUE)
Mset[,"varpos"] <- Mset[,indexname1]
Mset[,"varneg"] <- Mset[,indexname1]
Mset$varpos[which(Mset$varpos > 0)] <- rep(0, length(which(Mset$varpos > 0)))
Mset$varneg[which(Mset$varneg < 0)] <- rep(0, length(which(Mset$varneg < 0)))
commonxM <- Mset$time %in% indices$time
commonxI <- indices$time %in% Mset$time
subsetM <- Mset[commonxM,indexname1]
subsetI <- indices[commonxI,indexname2]
iok1 <- !is.na(subsetM) & !is.na(subsetI)
r=cor(subsetM[iok1,],subsetI[iok1,])
legname <- paste(indexname2," (r = ",toString(floor(r*100)/100),")",sep="")
p <- ggplot() +
  geom_area(data = Mset, aes(x=time, y=varpos), fill="blue") +
  geom_area(data = Mset, aes(x=time, y=varneg), fill="red") +
  geom_line(data = indices, aes(x=time, y=eval(as.name(indexname2))/(max2/max1), color=legname),size=0.5) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*(max2/max1),name = indexname2)) +
  scale_x_date(limits = c(min(Mset$time),max(Mset$time)),expand=c(0,0)) +
  labs(x  ="Time", y = indexname1) +
  theme(legend.position=c(0.15, 0.95),legend.key.width=unit(3,"line"),legend.background = element_rect(fill="transparent"),
        legend.text=element_text(size=11)) +
  scale_colour_manual("", 
                      breaks = c(legname),
                      values = c("black"))
p2=plotly_build(p)   
style(p2) %>% layout( legend = list(x = 0.01, y = 1) )   # needed because ggplotly removes the legend position